!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!!
!! Collection of Runge-Kutta time integrators
!!
!! \n
!!
!! Provide various RK integrators.
!<----------------------------------------------------------------------
module mod_rk

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_time_integrators_base, only: &
   mod_time_integrators_base_initialized, &
   c_ode

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_rk_constructor, &
   mod_rk_destructor,  &
   mod_rk_initialized, &
   ee_init,      ee_step,      ee_clean,      &
   hm_init,      hm_step,      hm_clean,      &
   ssprk54_init, ssprk54_step, ssprk54_clean

 private

!-----------------------------------------------------------------------

! Module variables

 ! public members
 logical, protected ::               &
   mod_rk_initialized = .false.

 ! private members

 !> time step
 real(wp) :: dt

 character(len=*), parameter :: &
   this_mod_name = 'mod_rk'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_rk_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized             .eqv..false.) .or. &
       (mod_kinds_initialized                .eqv..false.) .or. &
       (mod_time_integrators_base_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_rk_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_rk_initialized = .true.
 end subroutine mod_rk_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_rk_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_rk_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_rk_initialized = .false.
 end subroutine mod_rk_destructor

!-----------------------------------------------------------------------
 
 !> Explicit Euler initialization
 subroutine ee_init(dt_step)
  real(wp), intent(in) :: dt_step
 
  character(len=*), parameter :: &
    this_sub_name = 'ee_init'

   dt = dt_step
 
 end subroutine ee_init
 
!-----------------------------------------------------------------------
 
 !> Explicit Euler step
 subroutine ee_step(unew,t,unow,tmp)
  real(wp),     intent(in)    :: t
  class(c_ode), intent(in)    :: unow
  class(c_ode), intent(inout) :: unew
  !> \todo this argument should become an allocatable module variable
  class(c_ode), intent(inout) :: tmp(:)
 
  character(len=*), parameter :: &
    this_sub_name = 'ee_step'

   ! compute the tendency
   call tmp(1)%rhs(t,unow)
   ! compute the time step
   ! unew = unow + dt*tens
   call unew%op(unow,dt,tmp(1),tmp(2))
 
 end subroutine ee_step
 
!-----------------------------------------------------------------------
 
 !> Explicit Euler clean-up
 pure subroutine ee_clean()
  character(len=*), parameter :: &
    this_sub_name = 'ee_clean'

  ! nothing to do
 
 end subroutine ee_clean
 
!-----------------------------------------------------------------------

 !> Heun method initialization
 subroutine hm_init(dt_step)
  real(wp), intent(in) :: dt_step
 
  character(len=*), parameter :: &
    this_sub_name = 'hm_init'

   dt = dt_step
 
 end subroutine hm_init
 
!-----------------------------------------------------------------------
 
 !> Heun method step
 !!
 !! \f{displaymath}{
 !!  \begin{array}{rcl}
 !!   \displaystyle
 !!   u^{(1)} & = & \displaystyle u^n + \Delta t\,f(u^n) \\[3mm]
 !!   \displaystyle
 !!   u^{n+1} & = & \displaystyle u^n + \frac{\Delta t}{2} \left[
 !!   f(u^n) + f(u^{(1)}) \right].
 !!  \end{array}
 !! \f}
 !<
 subroutine hm_step(unew,t,unow,tmp)
  real(wp),     intent(in)    :: t
  class(c_ode), intent(in)    :: unow
  class(c_ode), intent(inout) :: unew
  !> \todo this argument should become a module variable
  class(c_ode), intent(inout) :: tmp(:)
 
  character(len=*), parameter :: &
    this_sub_name = 'hm_step'

   ! compute the tendency
   call tmp(1)%rhs(t,unow)
   ! compute the first time step
   call tmp(3)%op(unow,dt,tmp(1),tmp(2))

   ! compute the tendency
   call tmp(2)%rhs(t+dt,tmp(3))
   ! compute the second time step
   call tmp(3)%add(tmp(1),tmp(2))
   call unew%op(unow,0.5_wp*dt,tmp(3),tmp(1))

 end subroutine hm_step
 
!-----------------------------------------------------------------------
 
 !> Heun method clean-up
 pure subroutine hm_clean()
  character(len=*), parameter :: &
    this_sub_name = 'hm_clean'

  ! nothing to do
 
 end subroutine hm_clean
 
!-----------------------------------------------------------------------

 !> SSPRK54 method initialization
 subroutine ssprk54_init(dt_step)
  real(wp), intent(in) :: dt_step
 
  character(len=*), parameter :: &
    this_sub_name = 'ssprk54_init'

   dt = dt_step
 
 end subroutine ssprk54_init
 
!-----------------------------------------------------------------------
 
 !> SSPRK54 method step
 !!
 !! This is the method discussed in Table A.2 of <a
 !! href="http://dx.doi.org/10.1137/S0036142901389025">[Spiteri,
 !! Ruuth, SIAM J. Num. An., 2002]</a> and <a
 !! href="http://dx.doi.org/10.1007/BF01933264">[Kraaijevanger, BIT
 !! Numerical Mathematics, 1991]</a>.
 !!
 !! In the implementation, we make sure that the consistency condition
 !! \f{displaymath}{
 !! \sum_{k=0}^{i-1} a_{ik} = 1
 !! \f} is verified up to the working precision. Moreover, we notice
 !! that the intermediate time levels can be deduced (as fractions of
 !! \f$\Delta t\f$), setting \f$u^n=0\f$ and \f$f = 1\f$, so that
 !! \f{displaymath}{
 !!  \begin{array}{rclcl}
 !!    t^{(1)} & = & b_{10} & \approx & 0.39175 \\
 !!    t^{(2)} & = & a_{21}\,t^{(1)}+b_{21} & \approx & 0.58608 \\
 !!    t^{(3)} & = & a_{32}\,t^{(2)}+b_{32} & \approx & 0.47454 \\
 !!    t^{(4)} & = & a_{43}\,t^{(3)}+b_{43} & \approx & 0.93501 \\
 !!    t^{(5)} & = & a_{52}\,t^{(2)} + a_{53}\,t^{(3)} +
 !!    a_{54}\,t^{(4)} + b_{53}+b_{54} & = & 1.
 !!  \end{array}
 !! \f}
 !! Notice that \f$b_{54}\f$ is in fact computed using the constraint
 !! \f$t^{(5)} = 1\f$.
 !<
 subroutine ssprk54_step(unew,t,unow,tmp)
  real(wp),     intent(in)    :: t
  class(c_ode), intent(in)    :: unow
  class(c_ode), intent(inout) :: unew
  !> \todo this argument should become a module variable
  class(c_ode), intent(inout) :: tmp(:)

  real(wp), parameter :: &
    one = 1.0_wp
  real(wp), parameter :: &
    a10 = one,                                             &
    a20 = 0.44437049406734_wp, a21 = one-a20,              &
    a30 = 0.62010185138540_wp, a32 = one-a30,              &
    a40 = 0.17807995410773_wp, a43 = one-a40,              &
    a50 = 0.00683325884039_wp, a52 = 0.51723167208978_wp,  &
          a53 = 0.12759831133288_wp, a54 = one-(a50+a52+a53)
  real(wp), parameter :: &
    b10 = 0.39175222700392_wp, &
    b21 = 0.36841059262959_wp, &
    b32 = 0.25189177424738_wp, &
    b43 = 0.54497475021237_wp, &
    b53 = 0.08460416338212_wp ! b54 defined later, after the ta*
  real(wp), parameter :: &
    ta1 = b10,           &
    ta2 = a21*ta1 + b21, &
    ta3 = a32*ta2 + b32, &
    ta4 = a43*ta3 + b43, &
    ta5 = one
  real(wp), parameter :: &
    b54 = one-(a52*ta2+a53*ta3+a54*ta4+b53)

  character(len=*), parameter :: &
    this_sub_name = 'ssprk54_step'

   ! we use the temporaries as follows:
   !  1:4 stages (stage 5 is stored in unew)
   !  5:6 tendencies
   !  7:8 work array

   ! stage 1
   call tmp(5)%rhs(t,unow)
   call tmp(7)%stimes(a10,unow)
   call tmp(1)%op(tmp(7),dt*b10,tmp(5),unew)

   ! stage 2
   call tmp(5)%rhs(t+ta1*dt,tmp(1))
   call tmp(7)%stimes(a20,unow)
   call tmp(8)%op(tmp(7),  a21 ,tmp(1),unew)
   call tmp(2)%op(tmp(8),dt*b21,tmp(5),unew)

   ! stage 3
   call tmp(5)%rhs(t+ta2*dt,tmp(2))
   call tmp(7)%stimes(a30,unow)
   call tmp(8)%op(tmp(7),  a32 ,tmp(2),unew)
   call tmp(3)%op(tmp(8),dt*b32,tmp(5),unew)

   ! stage 4
   call tmp(5)%rhs(t+ta3*dt,tmp(3))
   call tmp(7)%stimes(a40,unow)
   call tmp(8)%op(tmp(7),  a43 ,tmp(3),unew)
   call tmp(4)%op(tmp(8),dt*b43,tmp(5),unew)

   ! stage 5
   call tmp(6)%rhs(t+ta4*dt,tmp(4))
   call tmp(7)%stimes(a50,unow)
   call tmp(8)%op(tmp(7),  a52 ,tmp(2),unew)
   call tmp(7)%op(tmp(8),  a53 ,tmp(3),unew)
   call tmp(8)%op(tmp(7),  a54 ,tmp(4),unew)
   call tmp(7)%op(tmp(8),dt*b53,tmp(5),unew)
   call   unew%op(tmp(7),dt*b54,tmp(6),tmp(1))

 end subroutine ssprk54_step
 
!-----------------------------------------------------------------------
 
 !> SSPRK54 method celan-up
 pure subroutine ssprk54_clean()
  character(len=*), parameter :: &
    this_sub_name = 'ssprk54_clean'

  ! nothing to do
 
 end subroutine ssprk54_clean
 
!-----------------------------------------------------------------------

end module mod_rk

